Reducing reflections from mesh refinement interfaces in numerical relativity 
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Full interpretation of data from gravitational wave observations will require accurate numerical 
simulations of source systems, particularly binary black hole mergers. A leading approach to improv- 
ing accuracy in numerical relativity simulations of black hole systems is through fixed or adaptive 
mesh refinement techniques. We describe a manifestation of numerical interface truncation error 
which appears as slowly converging, artificial reflections from refinement boundaries in a broad class 
of mesh refinement implementations, potentially compromising the effectiveness of mesh refinement 
techniques for some numerical relativity applications if left untreated. We elucidate this numerical 
effect by presenting a model problem which exhibits the phenomenon, but which is simple enough 
that its numerical error can be understood analytically. Our analysis shows that the effect is caused 
by variations in finite differencing error generated across low and high resolution regions, and that its 
slow convergence is caused by the presence of dramatic speed differences among propagation modes 
typical of 3+1 relativity. Lastly, we resolve the problem, presenting a class of finite differencing 
stencil modifications, termed mesh-adapted differencing (MAD), which eliminate this pathology in 
both our model problem and in numerical relativity examples. 

PACS numbers: 02.60. Jh,02.70.Bf,04.25.Dm,04.30.Db,04.70.-s 



I. INTRODUCTION 

Recent years have seen a dramatic rise in opportunities 
for observing strong-field gravitational dynamics. New 
observations of dense black-hole-like objects, at stellar, 
intermediate, and supermassive scales are increasingly 
frequent. Anticipated gravitational wave observations by 
ground-based and space-based detectors are expected to 
capture information about these objects at moments of 
the strongest gravitational interactions [ij, |2( . Interpre- 
tation of data from any such observations will depend 
on theoretical modeling of the strong-field interactions of 
dense black-hole-like objects in the process of generating 
gravitational radiation. General Relativity is the stan- 
dard model for describing gravitational interactions and 
wave generation. However, the predictions of General 
Relativity for such cases are not yet fully understood, 
and will depend on 3-D numerical relativity computer 
simulations Q- 

While numerical relativity has progressed markedly in 
recent years 0|, significant improvements in the fidelity 
of models for events such as binary black hole coalescence 
will be essential for the full interpretation of upcoming 
observations. There are many facets to the problem of 
improving such simulation,including optimally formu- 
lating Einstein's equations 0, IE 0> El j properly handling- 
boundaries, handling black hole singularities, handling 
constraint violations, and making judicious gauge choices 
. There are also basic numerical issues concerning how 
to, with finite resources, perform such high-fidelity 3-D 
simulations with strong short-wavelength gravitational 
features near the sources, and weak but critical long- 
wavelength gravitational (radiation) features emerging in 
a large domain. Approaches to this latter class of issues 
include, developing higher order accurate finite differ- 
encing methods El, spectral methods mesh refine- 
ment techniques |12lll3lll4l | and other forms of numerical 



patching techniques 0, El • 

We focus here on resolving a limitation which has 
arisen in our work on numerical relativity simulations 
of binary black hole systems with mesh refinement tech- 
niques, but which may have analogues in other numerical 
patching treatments as well. Mesh refinement techniques 
divide the computational domain into regions with sep- 
arate computational grids which can be of higher res- 
olution in some regions than others. Such approaches 
involve mesh-structure interfaces across which the de- 
tails of the finite differencing treatment suddenly change. 
Inevitably, these interfaces contribute to computational 
error, manifesting such effect as "reflections" off the in- 
terfaces. Considerable attention is given to implement- 
ing "clean" interfaces, which generate small error, com- 
pared to error generated in the bulk regions. At mini- 
mum, this requires interface-induced error to converge at 
least as rapidly as the bulk error. For some classes of 
black hole evolutions, with non- vanishing shift advection 
terms, we have typically observed large error propagating 
from mesh interfaces. Understanding and resolving this 
problem forms the focus of this paper, and we demon- 
strate a solution, mesh-adapted differencing (MAD) , that 
works. 



II. INTERFACE PERFORMANCE 

The pertinent features of our numerical scheme are as 
follows. We are solving the 3+1 BSSN formulation of 
Einstein's equations [E ll2L UtI ITsI ] . Our gauge condition 
is numerically determined, typically with some variation 
of the 1+log slicing and hyperbolic Gamma-driver shift 
evolution equations . We integrate in time with the it- 
erative Crank-Nicholson method [l9| . All spatial deriva- 
tives are computed by second order accurate, centered 
differencing, save for advection derivatives, for which we 
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FIG. 1: Convergence plot for the Hamiltonian constraint error 
H at time t = 4M. There is a single puncture black hole cen- 
tered at the origin, and refinement boundaries at \x%\ = 1M, 
2M, and 4M. The finest grid spacing hf for each simulation 
is indicated in the figure, and the higher resolution is multi- 
plied by a factor of 4. For errors demonstrating second order 
convergence the curves should superpose. 



use second order accurate upwinded differencing. 

We are particularly interested in simulating gravita- 
tional radiation generated in black hole collisions. To 
resolve the black hole sources adequately, whilst pushing 
the computational grid boundary sufficiently far away, 
we use fixed mesh refinement (FMR) , as implemented by 
a software package for this purpose called PARAMESH 
[20| . With this implementation, the resolutions of two ad- 
jacent refinement regions always differ by a factor of two; 
i.e. the grid-spacing of the coarser region is double what 
it is in the finer region. "Ghostzones" or "guardcells" , 
typically two layers, are required to provide buffering be- 
tween refinement levels. These guardcells are filled in by 
interpolation. 

Whether the simulation performs adequately in the 
presence of refinement interfaces is a question of particu- 
lar concern to us. Inevitably, refinement boundaries are 
sources of numerical errors, although these reflections are 
often satisfactorily convergent and negligible. An excep- 
tion has plagued us in the case of a non-negligible shift, 
/3 1 , at refinement boundaries. In this case we observe a 
reflection pulse that propagates at a velocity of — /3 1 . 

An example can be seen in the case of single 
Schwarzschild black hole (in isotropic Schwarzschild co- 
ordinates) centered in a nested-box arrangement of re- 
finement regions. The Hamiltonian constraint provides a 
measure of the error in the runs and is plotted in Fig. ^ 
These "reflected" error waves, or "bumps" , notably prop- 
agating toward the black hole (at x = 0) from each inter- 
face, are strongest when the refinement interface is near- 
est to the black hole (where (3 l is largest). We have also 
noticed that the reflections seem to originate coincidently 
with the passing of an initial gauge pulse through the in- 



terface. Such gauge pulses are typical in black hole sim- 
ulations with "1+log" type lapse conditions, and propa- 
gate in this case at \/2 times the speed of light asymptot- 
ically. In Fig^thc convergence of the bumps is indicated 
by comparing the error in the Hamiltonian constraint at 
two resolutions, a moderate resolution run of hf = M/32 
and run with uniformly doubled resolution, hf = Af/64, 
where hf refers to the resolution of the finest grid. The 
curves have been rescaled so that they should superpose 
if the errors are second order convergent. Disturbingly, 
the figure indicates that these errors not to converge at 
reasonable grid resolutions, M/32 and . Whether these 
reflection errors would converge manifestly at sufficiently 
high resolution is difficult to determine in our black hole 
applications of limited achievable grid size. 

In any case, the practical effect is that it is a challenge 
to effectively control these errors which may dominate the 
error around refinement boundaries in the strong field re- 
gion. Interestingly, however, we have not observed any 
adverse effect, such as poor convergence, imprinted by 
these bumps on such phenomena as gravitational radia- 
tion measured far from a binary black hole system |13| . 
However, other physical quantities of interest, such as the 
event horizon, seem likely to be more sensitive to these 
near-field errors. 



III. LINEARIZED BSSN 

To understand the source of the /3-speed error exhib- 
ited in the last section, we begin by considering a lin- 
earized BSSN system with 1+log slicing and hyperbolic 
Gamma-driver shift. The system of equations is: 
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where a = a — 1 , (3\ = /3 l — f3g, (with /3q assumed spatially 
uniform for simplicity), and hij = 7y — 8%j. 

For a problem that varies only in one dimension, with 
no initial transverse components, if we assume plane- 
wave solutions (generalizable by Fourier analysis), then 
the above system of equations can be written in the form: 



dAu) = ikU\u) 



(9) 



3 



where 



^i(kx~ujt) 



(10) 



with a, 

B x , fa, 



B, fa, 
d>, K, 



I a \ 

B 

fa 


K 

h 

A 

v f / 

-ftT, /i, A, and T the amplitudes of a, 
, j4 xx , and respectively, and 

















_2_ 

ik 








o \ 










Aik 
3 





4 
3 








-fa 







3 
4ifc 




























1 
6 


-fa 


1_ 

6ik 













—ik 











-fa 



















4 
3 








-fa 


__2_ 

ik 







2ik 
3 








Aik 
3 





ik 
2 


-fa 


2 

— 3 

-fa J 


V 








iik 
3 





4 
3 









The eigenvalues are 
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_3/?o + 
and 



— |A) — |VPo +4- If I A) is the eigenvector associ- 
ated with eigenvalue A, and (A| is defined such that 
Ex|A)(A| = I, then 



M = ^A|A)(A| 



(12) 



By substituting Eq. (|12f> into Eq. © the system of evolu- 
tion equations can be decomposed into a series of advec- 
tion terms, each associated with a characteristic velocity 
equal to one of the eigenvalues. Note that, assuming 
fa <C 1, the /3-speed mode is much slower than any of 
the other non-zero-speed modes. As we are particularly 
concerned with this mode, it is instructive to write the 
evolution equation thus: 
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+ A^|A)<A|u) 

This equation indicates that disturbances originating 
in 4> and can propagate in (f> and fty at /3-speed and, 
further, that this is the only means of /?-speed propaga- 
tion allowed by this system. Thus 4> is uniquely signifi- 
cant in both generating and propagating /3-speed modes. 
Note that Eq. i|13[l resembles the 4> evolution equation in 
the BSSN system. 
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FIG. 2: Convergence plot for the error in ip. There is a re- 
finement boundary at x = 50. The coarse grid spacing h of 
each simulation is indicated in the figure, and the medium 
and high resolution curves have been multiplied by 4 and 16 
respectively, as appropriate to demonstrate second order con- 
vergence. 



IV. A MODEL PROBLEM 

Motivated by the discussion of the last section, a simple 
model for the generation and propagation of the /3-speed 
modes is a one-dimensional advection problem with an 
additional driving term, 



tp(x, t) = fid x (p(x, t) + f(x, t), 



(14) 



In our numerical simulations, it appears that the re- 
flected "bumps" may be triggered by a rapidly propa- 
gating gauge pulse which propagates outward early in 
the simulations. For our model problem, which we will 
term "Bumpy" because of its most salient feature, we will 
drive the advection equation with a pulse propagating at 
speed v, significantly faster than the advection speed fa 
We thus let f{x,t) — f(x — vt) where v is larger than f3 
and both are assumed positive. This model equation is 
simple enough to be solved directly, 



(p(x,t) 



V + /3 



f(x+(t-t')P,t')dt' 

(F(x - vt) - F(x + fa)), (15) 



where d x F(x) — f(x). For definiteness we can take the 
driving term to be a Gaussian pulse, f(x) — exp(— x 2 ), 
implying F(x) — (-\/(7r)/2)erf(j;). For a localized pulse, 
such as this, the F(x + fa) term in the exact solution will 
be negligible in the region of interest, near x = xq. 

To test if this model problem is sufficient, we have 
numerically evolved Eq. Ijl4|l using a 1-D finite difference 
code with a resolution dx — h for x < xq and dx — 2h for 
x > xq, realizing a refinement boundary at x = xq = 50. 
We evolved over the domain < x < 100 with periodic 
boundary conditions, using three-point upwind finite dif- 
ference stencil and a mesh refinement scheme similar to 
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that in our numerical simulations with PARAMESH. In 
Fig. |21 we show the errors in our evolution variable ip in 
the vicinity of the refinement interface for runs with sev- 
eral resolutions, h — 5/64,5/128 and 5/256 . In each 
case there clearly a reflection error bump propagating to 
the left, away from the refinement interface, which we 
find propagates at speed p. As before the errors have 
been rescaled so that that second order converging error 
should superpose. Comparing the two lower resolution 
run, we do not see good superposition, and the peaks 
appear to converge at roughly first order. The compar- 
ative appearance of these errors is similar to that seen 
for black hole simulations in Fig.^ However, running at 
higher resolutions, which is readily allowed by this 1-D 
model, we find that the reflected error is indeed second 
order convergent, as is suggested by comparing the two 
higher resolution runs in Fig.0 This convergence is man- 
ifest only at relatively high resolution (high relative to the 
wavelength - and higher than can be easily achieved with 
respect to a gravitational wavelength in a binary black 
hole run). Further experimentation has indicated that 
these results are not strongly affected by variations in 
the time-integration method or by changing among mesh 
refinement interfacing schemes which are consistent with 
the overall second order finite differencing accuracy. 



V. A FINITE DIFFERENCE ANALYSIS 

In this section we attempt to understand the numer- 
ical behavior of the Bumpy problem Eq. (|14fl by con- 
structing here an analytic model for the numerical error 
in our simulations. Noting that the time discretization, 
and the details of the interpolation scheme used in apply- 
ing refinement conditions seem not to be directly linked 
to the problematic error features we see, we model the 
error with as few assumptions as possible about these de- 
tails. We treat the numerical errors continuously in time, 
and we consider the effects of spatial finite differencing 
in terms of a continuous field ip(x,t, h) representing the 
numerical solution at (fine-grid) resolution dx — h. We 
then expand ip in orders of h, 



<p(x, t, h) = tp e (x, t) + h 2 cp 2 {x, t) + 0{h 3 ), 



(16) 



Where cp e (x,t) is understood to be the exact solution 
Eq. I|15|) ■ We consider the effect of our finite differencing 
scheme by replacing the spatial derivative appearing in 
Eq. i|14|) with a suitable finite difference operator D^. 
Thus, 



tp(x, t, h) = /3D h tp(x, t, h) + f(x, t). 



(17) 



As above, we include a refinement jump in the grid at 
x = xq, represented here by applying a coarser version 
of the finite difference stencil in the x > xq part of the 
spatial domain. This will be consistent with any refine- 
ment interface algorithm which applies the same finite 
difference stencil in both the coarse and fine regions, and 



applies a interpolative guard-cell filling algorithm at the 
interfaces which leads to finite differences at the inter- 
face which are consistently second order accurate. On a 
uniform grid the second order error term of a finite first 
derivative is generally proportional to the third derivative 
of the field, 

D h <p(x, t, h) = d x (p(x, t, h) + e 2 h 2 d 3 x <p{x, t, h) + 0(h 3 ). 

For the specific upwind differencing operator used in Sec- 
tion ]IV\ the stencil of which is, 



Dh 



2Eh — 2^ 2h 
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where E is the spatial translation operator defined such 
that Eh f{x] — f(x + h), the constant error coefficient 
turns out to be e 2 = —1/3. Including the refinement 
jump, we have 

D h <p = D h ip + e(x-x )(D 2h -D h )<p (19) 
= d x <p + (1 + 30{x - x ))e 2 h 2 dlip + 0(h 3 ). 
Substituting into Eq. i|17|) . and rearranging, yields 

h 2 ip 2 (x,t,h) = [-cp e (x,t) + fid x ip e (x,t) +f(x,t)] 

+h 2 /3[d x ip 2 (x,t 1 h) + 

e 2 (l + 3Q{x-x ))d*ip e (x,t)] 

+0(h 3 ). (20) 

Then noting that the first term vanishes, and taking the 
limit h — > 0, we derive 

<pa(x,t) = I3d x ^ 2 {x,t)+ (21) 
[3e 2 (l + 30(x ~ x ))d 3 x ip e (x,t), 

where we have used the notation ip 2 (x,t) = ip 2 (x,t,0). 
This is our model for the generation and propagation of 
finite differencing error in the numerical model problem 
in Section Hvl 

Since Eq. I|21|l is of the same form as Eq. I|21|l we can 
solve it the like fashion. We leave out the negligible F(x+ 
(3t) in ip e , substituting 

f{x,t) = (3e 2 {l + 3O(x-x ))d 3 v e (x,t) 

~ ~rTRa + 3G(x-xoMF{x-vt) 
V + p 

into the first line of Eq. I|15(l . Then, performing the in- 
tegral with careful attention to the presence of the step 
function, we get the solution 

(p 2 (x,t) = (s(x - vt) - s(x + (3t))(l + 30{x - x )) 

+3s(^(x-x +f3(t-—)) x 
p v 

(Q(x - x + [3t) - Q(x - x )) (22) 
Pe 2 



where, 
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■ exp (— x 2 ). 
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FIG. 3: Numerical error in <p, immediately after a pulse inci- 
dent on the refinement interface at x — 50 has passed through, 
generating a pulse of transmitted error and a contracted pulse 
of reflected error. The error is second order convergent and in 
agreement with the analytic prediction of the numerical error, 
as indicated. 

As in Section llVl the s(x+j3t) is negligible in the relevant 
region near the refinement interface. Thus the first term 
in Eq. 1)22(1 is effectively the differencing error associated 
with the upsweep in (p e which propagates across the grid 
at speed v. This part grows four times larger in the coarse 
x > xo region. The second term is quite interesting, as it 
propagates in the reverse direction at speed (3. It has the 
same s(x) shape as the forward propagating component, 
but it is reversed, and contracted by a factor (3/v. It 
is timed to originate at the interface as the first pulse 
crosses. The two terms combine make the full solution 
continuous at the interface. Heuristically, one could say 
that the second term is caused by the discontinuity in 
the differencing error, generated in order to produce a 
regular solution, and then, necessarily advecting away as 
required by the original model equation Eq. (|14|) . It is 
contracted because it propagates at a different speed than 
the first term, but their time dependences must match at 
the interface. 

Note that although we concretely consider second or- 
der finite differencing here, the finite differencing order is 
largely irrelevant in this analysis. The calculation can be 
directly adapted for the leading order in a higher order 
finite differencing by changing a few coefficients. 

We check that this analysis accurately describes the 
numerical errors in our Bumpy model evolutions by com- 
paring the predicted error with the numerical error re- 
sults. As Fig. 13 shows, the prediction of Eq. (|22|l agrees 
to high precision with the numerical simulation errors 
realized in high resolution runs. The figure shows both 
the rightward propagating wave which travels at velocity 
v with the driving pulse, and the shorter-wavelength re- 
flected wave. The reason for the slow convergence in our 
Bumpy simulations is now clear. Because the reflection 
must propagate at a significantly slower speed than the 
pulse which generated it, while their frequencies must 
match, the reflected pulse necessarily has a shortened 



wavelength. The key ingredients in this producing this 
effect are that we are simulating a system with strongly 
mismatched propagation speeds, and with the potential 
for generating mode-mixing reflections off numerical fea- 
tures such as our refinement interfaces such. In such a 
situation, with resolution just sufficient to accurately re- 
solve the long-wavelength features in the solution, the 
shortened wavelength of the numerical reflections makes 
them unresolvable without a significant increase in the 
resolution (say by a factor oiv/f3). We thus have reason 
to expect that such an effect is indeed the source of the 
bumpy reflection in our black hole runs. We will verify 
this by attempting to eliminating this type of error. 

VI. IMPROVING THE DIFFERENCING 
STENCILS 

Slowly converging errors of the type exposed in Sec- 
tions now seem likely to occur under fairly general 
circumstances. This explains our difficulties in many 
attempts to eliminate these effects by tweaking the in- 
terface conditions in various ways. Our analysis sug- 
gests that the only ways to avoid these problems, aside 
from eliminating the slowly propagating modes by set- 
ting (3 l = 0, is to remove the mode-mixing reflections at 
interfaces. The results of the last section clearly show 
that the reflected error in the Bumpy system is over- 
whelmingly due to the discontinuity in the differencing 
operator. More specifically, the reflection is related to 
the discontinuity in the second order truncation error. 
This observation suggests a solution. 

Consider using modified differencing stencils such that 
the coefficient of the second order truncation error in the 
fine grid region is multiplied by a constant factor qo while 
the coefficient of the second order truncation error in the 
coarse grid region is multiplied by a constant factor q\. 
Then the previously constant coefficient of the last 
section becomes: 

62 -> [qo + [qi - qo)Q(x - x )}e 2 (23) 

Repeating the steps of the last section for this new 
truncation error, the solution for the second order error 
in if becomes, 

tpa(x,t) = (s(x - vt) - s(x + (3t))(q + (Aqi - q )Q{x - x )) 

+(M - q )s(-^(x -x + 0(t - — )) x 

p v 

(Q(x - x + fit) - Q(x - x )) (24) 

Thus for the spatially blueshifted, reflected error, we now 
obtain, 

Vref = (Ml - qo)s(~(x ~ -T + (3{t - ^))) (25) 
p V 

One possibility would be to eliminate the leading or- 
der reflection error simply by choosing qo = q± = 0. 
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This choice corresponds to using higher order stencils, 
say third or fourth order accurate, throughout the entire 
grid. Higher order differencing methods are clearly valu- 
able in numerical relativity simulations |fOf. However, 
we will focus here on identifying a minimal way to real- 
ize effective second order convergence in a second order 
convergent finite differencing scheme. Our result should 
generalize to arbitrary order. 

Note that with the choice qo = 1 and q\ = \, the 
reflected error vanishes. More generally, of course, any 
choice of qo and qi that makes the truncation error con- 
tinuous across the refinement boundary will remove the 
second order reflection. In particular, the choice 
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ho 



(26) 



will work in the nth refinement region of a grid with an 
arbitrary number of refinement levels, where ho is as- 
sumed to be the grid-spacing in the finest region. Wc 
call such a mesh-adapted differencing scheme MAD. 

A second order MAD stencil can be obtained simply by 
linearly combining a second order accurate stencil with 
a higher order accurate stencil, as follows: 



Dh n 



(i 



(27) 



where the superscripted numbers in square brackets rep- 
resent the order of accuracy of the differencing operator, 
and j > 2. (Of course, this stratagem can be readily gen- 
eralized for higher order MAD operators as well.) In the 
particular case of a second order upwinded stencil com- 
bined with a third order lopsided stencil, the resulting 
stencil has the form: 
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(28) 



We have implemented this MAD stencil in the numer- 
ical simulations of the Bumpy system. Fig. 0] shows the 
result for a run at a moderate resolution. With the MAD 
stencil the dominant reflected error is nearly eliminated. 



VII. RESULTS FOR BLACK HOLE 
EVOLUTIONS 

We can demonstrate that our analysis does indeed ex- 
plain the problematic reflections in our black hole sim- 
ulations, by verifying that the same remedies applied to 
fix our Bumpy model problem, i.e. higher order or MAD 
stenciling, also fix the problem in our black hole simula- 
tions. 

As in [13, we have found that third or higher order ac- 
curate differencing of advection terms is unstable in the 
vicinity of black hole punctures, if the stencil has any 



FIG. 4: Comparison of the error in tp reflected from the re- 
finement boundary at x — 50 in the case of a second order 
conventional stencil and a second order MAD stencil. 



points on the "downwind" side. If the 3rd or higher or- 
der accurate differencing stencil does not have any points 
on the downwind side, it requires three or more layers of 
guardcells to accommodate points on the upwind side, 
which is expensive memory- wise. Thus, although higher 
order spatial differencing should reduce reflection from 
refinement boundaries, it brings with it a new set of prob- 
lems to solve. 

Use of a second order accurate MAD stencil, as in 
Eq. I)28[l . for advection, avoids the above difficulties. As 
we generally locate the punctures within the finest grid 
regions, and the MAD stencil automatically reverts to 
conventional second order upwinding in this region, the 
advection derivative does not take any points on the 
downwind side in the vicinity of the puncture and we 
find that stability is maintained. 

In our Einstein solver, we have implemented a second 
order accurate MAD stencil for advection. For all non- 
advection derivatives, the MAD stencil is constructed 
from a linear combination of second order centered and 
fourth order centered stencils, as in Eq. (|27|l . As a result, 
reflection from refinement boundaries has been dramati- 
cally reduced in the case of a single puncture. Fig.[S]com- 
pares the Hamiltonian constraint error from a run with 
conventional stencils with the Hamiltonian constraint er- 
ror in runs with MAD stencils, demonstrating clear im- 
provement. The plot shows that the reflection error for 
our hf — M/32 case is much reduced with the MAD sten- 
cil. By comparing with the hf = M/64, MAD-stencil 
run, we also see that the remaining bump now converges 
away at second order or better. 

For simplicity we only show plots from single black hole 
simulations here, though we have applied MAD in more 
interesting binary black hole cases as well. We find that 
the improvements generalize naturally to these cases. In 
particular we find that errors in the Weyl scalars near 
refinement boundaries are reduced by at least an order 
of magnitude. 



7 




x(M) 

FIG. 5: Comparison of the Hamiltonian constraint H in the 
case of second order accurate, conventional differencing sten- 
cils and second order accurate, MAD stencils, as well as a 
demonstration of convergence in the latter case. There is a 
single puncture black hole centered at the origin, and refine- 
ment boundaries at \xA = 1M, 2M, and 4M. 



VIII. CONCLUSIONS 

We have studied a problem which occurs in numerical 
relativity simulations with nonvanishing shift with mesh 
refinement. These involve slow advection modes across a 
mesh interface, and produce slowly converging numerical 



reflections that propagate at the speed (3 of the advection. 
We have successfully modeled the problem analytically. 
Our analysis suggests that the the effect is a general con- 
sequence of discontinuous jumps in the finite differencing 
stencil in a problem with propagation modes of widely 
differing speeds. Our proposed solution, to adjust the fi- 
nite differencing stencils to make the leading order differ- 
encing error continuous across mesh interfaces, is shown 
to be effective in black hole simulations with mesh refine- 
ment, but may have wider application in numerical rel- 
ativity. In addition to grids of non-uniform refinement, 
mesh-adapted differencing may also be appropriate for 
grids with multiple coordinate patches, where discontin- 
uous differencing error can also be expected. 



Acknowledgments 

We are happy to thank J. David Brown for help- 
ful discussion. We gratefully acknowledge CPU time 
grants from the Commodity Cluster Computing Project 
(NASA-GSFC) and Project Columbia (NASA Advanced 
Supercomputing Division, NASA- Ames). This work was 
supported in part by NASA Space Sciences grant ATP02- 
0043-0056. JvM was also supported in part by the Re- 
search Associateship Programs Office of the National Re- 
search Council. 



[1] G. Gonzalez (LIGO Scientific), Pramana 63, 663 (2004). 
[2] K. Danzmann and A. Rudiger, Class. Quant. Grav. 20 
(2003). 

[3] B. Schutz, in The Astrophysics of Gravitational Wave 
Sources, edited by J. M. Centrella (AIP, Melville, NY, 
2003). 

[4] M. Alcubierre (2004), gr-qc/0412019. 

[5] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 

024007 (1999), gr-qc/9810065. 
[6] H. Friedrich, Class. Quantum Grav. 13, 1451 (1996). 
[7] L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, Phys. 

Rev. D 64, 064017 (2001), gr-qc/0105031. 
[8] C. Bona, T. Ledvinka, C. Palenzuela, and M. Zacek, 

Phys. Rev. D67, 104005 (2003), gr-qc/0302083. 
[9] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 

D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 

67, 084023 (2003). 
[10] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O. 

Lousto (2005), gr-qc/0505055. 
[11] L. E. Kidder, M. A. Scheel, S. A. Teukolsky, E. D. Carl- 
son, and G. B. Cook, Phys. Rev. D 62, 084032 (2000), 



gr-qc/0005056. 

[12] B. Imbiriba, J. Baker, D.-I. Choi, J. Centrella, D. R. 
Fiske, J. D. Brown, J. van Meter, and K. Olson, Phys. 
Rev. D 70, 124025 (2004), gr-qc/0403048. 

[13] D. R. Fiske, J. G. Baker, J. R. van Meter, D.-I. Choi, and 
J. M. Centrella, Phys. Rev. D (2005), gr-qc/0503100. 

[14] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quan- 
tum Grav. 21, 1465 (2004), gr-qc/0310042. 

[15] G. Calabrese and D. Neilsen (2004), gr-qc/0412109. 

[16] J. Thornburg, Class. Quant. Grav. 21, 3665 (2004), gr- 
qc/0404059. 

[17] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 
(1995). 

[18] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. 

Phys. Suppl. 90, 1 (1987). 
[19] S. Teukolsky, Phys. Rev. D 61, 087501 (2000), gr- 

qc/9909026. 

[20] P. MacNeice, K. M. Olson, C. Mobarry, R. deFainchtein, 
and C. Packer, Computer Physics Comm. 126, 330 
(2000). 



